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Abstract The numerical analysis of variational integrators relies on variational 
error analysis, which relates the order of accuracy of a variational integrator 
with the order of approximation of the exact discrete Lagrangian by a com- 
putable discrete Lagrangian. The exact discrete Lagrangian can either be char- 
acterized variationally, or in terms of Jacobi's solution of the Hamilton-Jacobi 
equation. These two characterizations lead to the Galerkin and shooting-based 
constructions for discrete Lagrangians, which depend on a choice of a numerical 
quadrature formula, together with either a finite-dimensional function space or 
a one-step method. We prove that the properties of the quadrature formula, 
finite-dimensional function space, and underlying one-step method determine 
the order of accuracy and momentum-conservation properties of the associated 
variational integrators. We also illustrate these systematic methods for con- 
structing variational integrators with numerical examples. 



1 Introduction 

Variational integrators are a particularly popular class of geometric numerical 
integrators, due to their symplectic and momentum conservation properties, 
and their good energy behavior. While most of the standard approaches for 
constructing variational integrators involve piecewise polynomial interpolants, 
and high-order quadrature formulas, we will show that existing techniques from 
approximation theory, numerical quadrature, and one-step methods, can be 
incorporated systematically into the variational integration framework. 

In this paper we will discuss two general approaches to constructing varia- 
tional integrators. The first of this is the Galerkin variational integrator con- 
struction, which relies on a variational characterization of the exact discrete 
Lagrangian, and involves the choice of a numerical quadrature formula, and a 
finite-dimensional function space. The second approach relies on the charac- 
terization of the exact discrete Lagrangian in terms of Jacobi's solution of the 
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Hamilton-Jacobi equation, and involves the choice of a numerical quadrature 
formula, and an underlying one-step method. 

1.1 Discrete Lagrangian Mechanics 

Discrete Lagrangian mechanics |21.] is based on a discrete analogue of Hamilton's 
principle, referred to as the discrete Hamilton's principle, 

SSd = 0, 

where the discrete action sum, Sd '■ Q"^^ R, is given by 

En— 1 
. Ld{q^,q^+l). 

The discrete Lagrangian, : Q x Q — > M, is a generating function of the 
symplectic flow, and is an approximation to the exact discrete Lagrangian, 

L^iqo,qi;h)^ / Liqoiit), qoiit))dt, (1) 
Jo 

where (701 (0) = qo, 901 (^) = 91 j Sind qoi satisfies the Euler-Lagrange equation in 
the time interval (0, ft.). The exact discrete Lagrangian is related to the Jacobi 
solution of the Hamilton-Jacobi equation. Alternatively, one can characterize 
the exact discrete Lagrangian in the following way, 

L^{qo,qi;h)= ext / L{q{t),q{t))dt. (2) 

9(o)=go,g('i)=9i 

As we will see in the remainder of the paper, these two characterizations of the 
exact discrete Lagrangian will lead to two general techniques for constructing 
variational integrators. 

The discrete variational principle then yields the discrete Eulei — Lagrange 
(DEL) equation, 

D2Ldiqk-i,qk) + DiLdiqk,qk+i) = 0, (3) 

which implicitly defines the discrete Lagrangian map F^^ : {qk-i,qk) '-^ 
{qk,qk+i) for initial conditions (go 7 91) that are sufficiently close to the diago- 
nal of Q X Q. This is equivalent to the implicit discrete Eulei — Lagrange 
(IDEL) equations, 

Pk = -DiLd{qk,qk+i), Pk+i = D2Ld{qk,qk+i), (4) 

which implicitly defines the discrete Hamiltonian map Fl^ : {qk,Pk) '-^ 
{qk+i,Pk+i), where the discrete Lagrangian is the Type I generating function 
of the symplectic transformation. 
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1.2 Desirable Properties of Variational Integrators 

We summarize the many geometric structure-preserving properties of varia- 
tional integrators, which account for their popularity in simulations involving 
mechanical systems. 

Symplecticity. Given a discrete Lagrangian L^, one obtains a discrete fiber 
derivative, ¥Ld : (90,91) '-^ (qo, —DiLd{qo,qi))- Variational integrators are 
symplectic, i.e, the pullback under FL^ of the canonical symplectic form 
n on the cotangent bundle T*Q, is preserved. Pushing-forward the dis- 
crete Euler-Lagrange equations yield a symplectic-partitioned Runge-Kutta 
method. 

Momentum Conservation. Noethcr's theorem states that if a Lagrangian is in- 
variant under the lifted action of a Lie group, the associated momentum is 
preserved by the flow. If a discrete Lagrangian is invariant under the diago- 
nal action of a symmetry group, a discrete Noether's theorem holds, and the 
discrete flow preserves the discrete momentum map. For PDEs with a uni- 
form spatial discretization, a backward error analysis implies approximate 
spatial momentum conservation [22]. 

Approximate Energy Conservation. While variational integrators do not exactly 
preserve energy, backward error analysis [l|; @; 0; [13] shows that it preserves 
a modified Hamiltonian that is close to the original Hamiltonian for expo- 
nentially long times. In practice, the energy error is bounded and does not 
exhibit a drift. This is the temporal analogue of the approximate momentum 
conservation result for PDEs, as energy is the momentum map associated 
with time invariance. 

Applicable to a Large Class of Problems. The discrete variational approach is 
very general, and allows for the construction of geometric structure-preserving 
numerical integrators for PDEs 18], nonsmooth collisions j5j, stochastic sys- 
tems 0, nonholonomic systems [3], and constrained systems Further- 
more, Dirac structures and mechanics allows for interconnections between 
Lagrangian systems, thereby providing a unified simulation framework for 
multiphysics systems. 

Generates a Large Class of Methods. A variational integrator can be constructed 
by choosing a finite-dimensional function space, and a numerical quadrature 
method [15| . By leveraging techniques from approximation theory, numer- 
ical analysis, and finite elements, one can construct variational integrators 
that are appropriate for problems that evolve on Lie grou ps j l2c il3.] and 
homogeneous spaces [3], or exhibit multiple timescales [isl: l25l|. 



1.3 Variational Error Analysis and Discrete Noether's Theorem 

The variational integrator approach to constructing symplectic integrators has 
a few important advantages from the point of view of numerical analysis. In 
particular, the task of proving properties of the discrete Lagrangian map F^^ : 
QxQ ^ QxQ reduces to verifying certain properties of the discrete Lagran gian 
instead. Here, we summarize the results from Theorems 1.3.3 and 2.3.1 of 21 1 
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that relate to the order of accuracy and momentum conservation properties of 
the variational integrator. 

Discrete N aether's theorem. Given a discrete Lagrangian Ld : Q y- Q ^ 
which is invariant under the diagonal action of a Lie group G on Q x Q, then 
the discrete Lagrangian momentum map, Jl^ : Q x Q — > 0*, given by 

JLMk,qk+i) ■ ^ = {-DiLd{qk,qk+i,^Q{qk)) 
is invariant under the discrete Lagrangian map, i.e., Jl^ o Fl^ — Jl^- 

Variational error analysis. The natural setting for analyzing the order of ac- 
curacy of a variational integrator is the variational error analysis framework 
introduced in In particular. Theorem 2.3.1 of [2l[ states that if a discrete 
Lagrangian, : Q x Q R, approximates the exact discrete Lagrangian, 
if : Q X Q ^ M, given in ^ and ^ to order p, i.e., 

Ld{qo, qi; h) = if (go, gi; h) + 

then the discrete Hamiltonian map, F^^ : {qk,Pk) '-^ iqk+i,Pk+i), viewed as a 
one-step method, is order p accurate. 



2 Galerkin Variational Integrators 

The variational characterization of the exact discrete Lagrangian ([2]), 

Ld{qo,qi;h) = ext / L{q{t),q{t))dt, 
qecH[o,h]:Q) Jo 

Q{0)=qO:q{h)=qi 

naturally leads to the Galerkin approach to constructing variational integrators. 
This involves replacing the infinite-dimensional function space C^([0, h], Q) with 
a suitably chosen finite-dimensional subspace, and the integral with a numerical 
quadrature formula. In particular, it is common to choose an interpolatory 
approximation space, so as to satisfy the boundary conditions at time and h 
in a straightforward manner. 

We now recall the construction of higher-order Galerkin variational integra- 
tors, as originally described in 21,]. Given a Lie group G, the associated state 
space is given by the tangent bundle TG. In addition, the dynamics on G is 
described by a Lagrangian, L : TG M. Given a time interval [0,h], the 
action map, & : C^([0, h],G) M, is given by 

&{g)= [ L{g{t),g{t))dt. 
Jo 

We approximate the action map, by numerical quadrature, to yield & : 

C^{[0,h],G) ^R, 



e^(ff) = /iV' bM9{c^h),g{cM, 

^ 
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where q G [0, 1], i = 1, . . . , s are the quadrature pomts, and bi are the quadra- 
ture weights. 

RecaU that the discrete Lagrangian should be an approximation of the form 
Ld{9Q,9i,h) Ki ext 6(5). 

geC^{[0,h],G), 

g{o)=go,gih)=gi 

If we restrict the extremization procedure to the subspace spanned by the in- 
terpolatory function that is parameterized by s -I- 1 internal points, ip : G"*^^ — >■ 
C^([0, /i], G), we obtain the following discrete Lagrangian, 

Ld{go,9i)= ext &{(p{g'';-)) 

g 

g°=go\g'=gi 

biL(Tip(g'^; Cih)). 

g°=go\g'=gi 



2.1 Galerkin Lie Group Variational Integrators 

A particularly novel and nontrivial application of this approach is the construc- 
tion of Lie group variational integrators. This is based on the underlying idea of 
Lie group integrators 8] , which is to express the update map of the numerical 
scheme in terms of the exponential map, 

91 = ffoexp(^oi) , 

and thereby reduce the problem to finding an appropriate Lie algebra element 
Coi G 0j such that the update scheme has the desired order of accuracy. This 
is a desirable reduction, as the Lie algebra is a vector space, and as such the 
interpolation of elements can be easily defined. In our construction, the inter- 
polatory method we use on the Lie group relies on interpolation at the level of 
the Lie algebra. 

Since the momentum conservation properties of variational integrators de- 
pend on the discrete Lagrangian being invariant under the diagonal action of 
the symmetry group, we will introduce a construction that yields a G-invariant 
discrete Lagrangian whenever the continuous Lagrangian is G-invariant. This 
is achieved through the use of G-equivariant interpolatory functions, and in 
particular, natural charts on G, which we will now discuss. This then leads to 
the issue of reduction of higher-order Lie group integrators. 

Group- equivariant interpolants. The interpolatory function is G-equivariant if 
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Proposition 1 If the interpolatory function ip{g'^]t) is G-equivariant, and the 
Lagrangian, L : TG M, is G-invariant, then the discrete Lagrangian, Ld '■ 
G X G — )■ M, given by 



Ldigo,gi)= ext hS^ biL{T(f{g'';Cih)), 

9° =9a\g' =a\ 

is G-invariant. 
Proof 

. b,L{Tip{g'';Cih)), 

9°=g9o;a°=ggi 

Es 
. b,L{TLp{gg'';Cih)), 

gg°=ggo;99°=ggi 



g =so;s =91 

h 



ext biL[Tip{g^ ]Cih)), 

9"eG; ' 



g =go;g =gi 
= Ld{go,gi), 

where we used the G-equivariance of the interpolatory function in the third 
equality, and the G-invariance of the Lagrangian in the forth equality. 

Remark 1 While G-equivariant interpolatory functions provide a computation- 
ally efficient method of constructing G-invariant discrete Lagrangians, we can 
construct a G-invariant discrete Lagrangian (when G is compact) by averag- 
ing an arbitrary discrete Lagrangian. In particular, given a discrete Lagrangian 
Ld : Q X Q ^ R, the averaged discrete Lagrangian, given by 



Ld{qo,qi) = -jy^ I Ld{gqo,gqi)dg 
1^1 JqeG 



geG 

is G-equivariant. Therefore, in the case of compact symmetry groups, a G- 
invariant discrete Lagrangian always exists. 

Natural charts. Following the construction in , we use the group exponential 
map at the identity, exp^ : g — >■ G, to construct a G-equivariant interpolatory 
function, and a higher-order discrete Lagrangian. As shown in Lemma [U this 
construction yields a G-invariant discrete Lagrangian if the Lagrangian itself is 
G-invariant. 

In a finite-dimensional Lie group G, expg is a local diffeomorphism, and thus 
there is an open neighborhood U C G of e such that exp^^ : U — > u C g. When 
the group acts on the left, we obtain a chart tpg : LgU u at g £ G by 



V'g — expg ^ oL 
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We would like to construct an interpolatory function that is described by a set 
of control points {g'^}t=Q in the group G at control times = do < di < ci2 < 
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. . . < ds-i < ds = 1. Our natural chart based at induces a set of control 
points ^'^ = "ipgoig") in the Lie algebra Q at the same control times. Let lu,s{t) 
denote the Lagrange polynomials associated with the control times d^, which 
yields an interpolating polynomial at the level of the Lie algebra, 

Applying ip^^ yields an interpolating curve in G of the form, 

ifig''; Th) = V-gO^ (X]l=o '^9°{9''%A'r)) , 

where (fi{d^h) = g" for z/ = 0, . . . , s. Furthermore, this interpolant is G-equivariant, 

as shown in the following Lemma. 

Proposition 2 The interpolatory function given by 

^{g"; Th) = ( ^l^^ i'go (gnlAr)) , 

is G-equivariant. 
Proof 

ifiigg^-^Th) = i'{glo)(^Yll=o^99'>(99''%A'r)) 

= Lggo expg (^Yll^o^^Pe\igg°)~\ggl%AT)) 

= LgLgO exp, [j2l=o^'^e\{9')-'9-'99nlAr)) 

= Largo'{Y.l=j9<9nlAr)) 

= Lgip{g";Th). 

Remark 2 In the proof that ip is G-equivariant, it was important that the base 
point for the chart should transform in the same way as the internal points 
g'\ As such, the interpolatory function will be G-cquivariant for a chart that it 
based at any one of the internal points g'^ that parameterize the function, but 
will not be G-equivariant if the chart is based at a fixed g G G. Without loss of 
generality, we will consider the case when the chart is based at the first point 

50- 

We will now consider a discrete Lagrangian based on the use of interpolation 
on a natural chart, which is given by 

^450,51) = ext 1 b^L{T^{{g•'}t=o■, Cih)) . 

a°=ao;g" =90^^91 

To further simplify the expression, we will express the extremal in terms of the 
Lie algebra elements ^'^ associated with the v-th control point. This relation is 
given by 

C = ^goi9n, 
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and the interpolated curve in the algebra is given by 
which is related to the curve in the group, 

The velocity ^ = g~^g is given by 

tirh) = g-'g{rh) = ^ Y^l^^ClAr)- 
Using the standard formula for the derivative of the exponential, 

exp = TeLexp(c) • dexpa^d^ , 

where 



dexp.=E„.o(n + l)!' 
we obtain the following expression for discrete Lagrangian, 

Ld{go,gi)= ext /lY]. ^i-^f^so exp(C(ci/i)), 
«°=0;l''=Vpo(si) 

Te^p{S{cih))Lgo • Teiexp(«(ci/»)) " dcxp^^^^^^^^ (^(Ci/l))) . 

More explicitly, we can compute the conditions on the Lie algebra elements for 
the expression above to be extremal. This implies that 

Ld{go,gi) = /iX]i=i ^'^{^90 exp(^(ci/i)), 



with ^° = 0, = tpgoigi), and the other Lie algebra elements implicitly defined 

by 

Es dL 
i=l ■^('^»'*)^exp(i(ci/i))igo • ^e-t'exp(€(ci/i)) " dexPg^j^^^^^j lv,s{Ci) 

1 dL ~ 

+ /^■^(CiM^exp(C(c./i))-^exp(4(c.ft)) " T'e -£'exp(C(c.?i)) " ddexPg^^^^,,;.) ^.^(Ci) 



for 1/ = 1, . . . , s — 1, and where 

n — 



n=o (n + 2)! ' 

This expression for the higher-order discrete Lagrangian, together with the 
discrete Euler-Lagrange equation, 

D2Ld{go,gi) + DiLd{gi,g2) = 0, 
yields a higher-order Lie group variational integrator. 
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2.2 Higher-Order Discrete Euler-Poincare Equations 

In this section, we will apply discrete Euler-Poincare reduction (see, for ex- 
ample, 0) to the Lie group variational integrator we derived previously, to 
construct a higher-order generalization of discrete Euler-Poincare reduction. 



Reduced discrete Lagrangian. We first proceed by computing an expression for 
the reduced discrete Lagrangian in the case when the Lagrangian is G- invariant. 
Recall that our discrete Lagrangian uses G-equivariant interpolation, which, 
when combined with the G-invariance of the Lagrangian, implies that the dis- 
crete Lagrangian is G-invariant as well. We compute the reduced discrete La- 
grangian, 

h{go^gi) = Ld{ga,g{) 

= Ld{e,g^^gi) 

= ext h'S^ biL(Leexp{^{cih)), 

C°=0;C = =log(3-isi) 



ext biLl expi^icih)), 



«°=0;« = =log(go"igi) 



T, 



eicxp(^(c.?i)) • dexp^d5(,^^)('?(cj/i)) 

Setting ^'^ = 0, and = \og{gQ^gi), we can solve the stationarity conditions 
for the other Lie algebra elements using the following implicit system 

of equations. 



dL 

— (c,ft.)reLexp(^(c,/i)) ■ dexpg^d^^^^^^ L,siCi) 

1 dL ~ 
+ 7^^(c»/i)7'eioxp(e(c.M) -ddexp^d^t,,,,, l.s{c,) 



where v = 1, . . . , s — 1. 

This expression for the reduced discrete Lagrangian is not fully satisfactory 
however, since it involves the Lagrangian, as opposed to the reduced Lagrangian. 
If we revisit the expression for the reduced discrete Lagrangian, 

kido^gi) = 

&d ^«-^(exp(C(ci^)),reLcxp(C(c./i))-dexp^d5(,^;,,(C(ciM)) > 

{O=0;e»=log'(g-igi) 

we find that by G-invariance of the Lagrangian, each of the terms in the sum- 
mation. 



Li exp(^(cift.)),reLexp(5(c./i)) ' dexp^j^^^ (^(q/i)) 



10 



Molvin Look, Tatiana Shingel 



can be replaced by 

;(dexp^d^^^^^^^(C(c,/i))) , 

where Z : g — M is the reduced Lagrangian given by 

Kv) = L{Lg-ig,TLg-ig) = L(e,ry), 

where rj = TLg-ig G g. 

From this observation, we have an expression for the reduced discrete La- 
grangian in terms of the reduced Lagrangian, 

Idigo^gi)^ ^ext /i^.^^ foj^dexp^j^^^^^^ (^(c,/i)) 

?°=0;C==log'(g-igi) 

As before, we set = 0, and = log{gQ^gi), and solve the stationarity con- 
ditions for the other Lie algebra elements {^''j^Z^ using the following implicit 
system of equations, 

= hY,^^^ b, Q-{c^h) ddexp,d,(.,,) l.siC'O 



where v = 1, . . . , s — 1. 



Discrete Euler-Poincare equations. As shown above, we have constructed a 
higher-order reduced discrete Lagrangian that depends on 

/fcfc+i = 9k9k+i- 

We will now recall the derivation of the discrete Euler-Poincare equations, 
introduced in 20]. The variations in fkk+i induced by variations in gk, gk+i are 
computed as follows, 

Sfkk+i = -g,:^Sgkgk-lgk+i +g,:^6gk+i 

= TRh^+A-9k^^9k + Mf^^_^, gk+iSgk+i) • 

Then, the variation in the discrete action sum is given by 



k=Q ^d{fkk+l)Sfkk+l 

^l^k=o '■'i(fkk+i)TRf,,+A-9k Sgk+ Mf,,^,gk+iSgk+i) 
= Er=i' [l'd{fk-ik)TRf,_,, Ad/,_,, -l'd{fkk+i)TRf,,^,] dk 



with variations of the form -dk — 9k^^9k- In computing the variation of the 
discrete action sum, we have collected terms involving the same variations, and 
used the fact that do = = 0. This yields the discrete Eulei — Poincare 
equation, 

Uh-ik)TRf,_,, Ad/,_,, ~Ufkk+i)TRf,,^, =0, fc = 1, . . . , TV - L 
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For ease of reference, we will recall the expressions from the previous discussion 
that define the higher-order reduced discrete Lagrangian, 

ld{fkk+l) = ^X]i=i ''''(^^^Pad^ccjfc) 



where 



and 



e = iog(/fcfe+i) , 

and the remaining Lie algebra elements {^''l^zi, are defined implicitly by 

for 1/ = 1, . . . , s — 1, and where 

ddexp^ = Er=o ■ 

When the discrete Euler Poincare equation is used in conjunction with the 
higher-order reduced discrete Lagrangian, we obtain the higher-order Euler- 
Poincare equations. 



^^{c,h) ddexp^d^^^^^j K^,{ci) 



2.3 Example: Lie Group Velocity Verlet 

We will now construct a Lie group analogue of the velocity Verlet method for 
the free rigid body. The velocity Verlet method can be derived from the context 
of discrete mechanics by considering the following discrete Lagrangian, 

Ld{qk,qk+i) = 2 

which corresponds to using a piecewise linear interpolant, and the trapezoidal 

rule to approximate the integral. 

In the case of the free rigid body, the Lagrangian is given by, 

L{R,R) = ]^njn^ = \^r{S{n)JiS{Qf\ . 

Here, Jd is a modified moment of inertia that is related to the usual moment of 

inertia by the relations, = 5(tr[J] /sxs — 2J), and J = tr[Jd] /sxs — Jd- From 
the kinematic relation Siji) = BFR, we have that, 

S{f2k) — R^Rk ~ Rk = Ti^k — hxs), 



qk, 



Qk+i — qk 
ll 



*+! - qk 



h 
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where Fk = Rj^Rk+i- Then, the discrete Lagrangian for the velocity Verlet 
method apphed to the free rigid body is given by, 

Ld{Rk:Rk+i) = 2 • ^^/^tr[(Ffc - h>c3fJd{Fk - /sxa)] 

- ■^tr[(/3x3 - Fk)Jd] , 

where in the second to last equality, we used the fact that tr[74_B] — tT[BA] , and 
in the last equality, we used the fact that F^. is an orthogonal matrix, is 
symmetric, and tr[A_B] — tr[i3-^^-^] . 

Recall that • 6R = —R^{5R)R^. Furthermore, the variation of Rk is 
given by, 

6Rk = Rkfik, 

where rjk G so (3) is a variation represented by a skew-symmetric matrix and 
vanishes aX k — Q and k = N . We may now compute the constrained variation 
of Fk = R]^Rk, which yields, 

6Fk = SR^Rk+i + RlSRk+i = rjkRlRk+i + RlRk+im+i = -VkFk + Fkrjk+i. 
Define the discrete action sum to be 

6d = Ld{Rk,Fk). 

Taking constrained variations of Fk yields, 

EN-l 1 
, „ T{ty[-r]k+iJdFk]+tr[7]kFkJd]}. 
k=a h 

Using the fact that the variations rjk vanish at the endpoints, we may reindex 
the sum to obtain, 

SBd = 22k=i Jl^i^lkiFkJd - JdFk-i)] . 

The discrete Hamilton's principle states that the variation of the discrete action 
sum should be zero for all variations that vanish at the endpoints. Since rjk is 
an arbitrary skew-symmetric matrix, for the discrete action sum to be zero, it 
is necessary for {FkJd — JdFk-i) to be symmetric, which is to say that, 

Fk+iJd — JdFk+i ~ JdFk + Fk -^d ~ 0. 

This implicit equation for Fk+\ in terms of Fk^ together with the reconstruction 
equation i?fc+i = RkFk, yields the Lie group analogue of the velocity Verlet 
method. 

In practice, in time marching the numerical solution, we need to solve the 
above equation for Fk+i G 50(3) given Fk- This equation is linear in Fk+i, but 
it is implicit due to the nonlinear constraint F/^^^Fk+i — /sxs- Since JdFk — 
Fk Jd is a skew-symmetric matrix, it may be represented as S{g), where g £ R^, 
which reduces the equation to the form, 

FJd ' JdF^ = S{g). (5) 

We now introduce two iterative approaches to solve ([5|) numerically. 
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Exponential map. An element of a Lie group can be expressed as the exponential 
of an element of its Lie algebra, so F £ SO (3) can be expressed as an exponential 
of S{f) G so (3) for some vector / £ R^. The exponential can be written in closed 
form, using Rodrigues' formula, 

F ^ exp Sif) = /3X3 + + '-^^Siff. (6) 

Substituting ([6|) into ([5]), we obtain 

Si9) - ^ W) + X J/). 

11/11 li/f 

Thus, ^ is converted into the equivalent vector equation g — G{f), where 
G : is given by 

1/11 ii/f 

We use the Newton method to solve g = G{f), which gives the iteration 

f^+l - /. + VG(/,)-' {g - G(/,)) . (7) 

We iterate until — G(/i)|| < e for a small tolerance e > 0. The Jacobian 
VG(/) in ([7|) can be expressed as 

ii/f 11/11 

sin 11/11 11/11 -2(1 -cos 11/11) (.^jr.rT 

11/11' 

1 — COS 



' ll^l|2 {-SiJf) + Sif)J}. 

Cayley transformation. Similarly, given /c G R'^, the Cayley transformation is 
a local diffeomorphism that maps S{fc) G so(3) to F e S0(3), where 

F = cay S{f,) - (/3X3 + S{fMhx3 ~ S{f,))-\ (8) 

Substituting ^ into (O, we obtain a vector equation Gdfc) — equivalent to 

m 

Gcifc) =9 + 9xfc + i9^fc)fc - 2Jf, = 0, (9) 

and its Jacobian VGc(/c) is written as 

VGc(/c) = Sig) + (ff^/c)/3x3 + fcg^ - 2J. 

Then, ([9]) is solved by using Newton's iteration (O, and the rotation matrix is 
obtained by the Cayley transformation. 

For both methods, numerical experiments show that 2 or 3 iterations are 
sufficient to achieve a tolerance of e = 10^^^. Numerical iteration with the 
Cayley transformation is a faster by a factor of 4-5 due to the simpler expressions 
in the iteration. It should be noted that since F = expS'(/) or F = cay S{fc), 
it is automatically a rotation matrix, even when the equation g = G{f) is not 
satisfied to machine precision. This computational approach is distinct from 
directly solving the implicit equation ([5]) with 9 variables and 6 constraints. 
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2.4 Numerical comparisons 



The Lie Group Velocity Verlet method is a second-order symplectic Lie group 
method, and it is natural to compare it to other second-order accurate methods 
which fail to preserve either the symplectic or Lie group structure: 

i. Explicit Midpoint Rule (RK): Preserves neither symplectic nor Lie group 
properties. 

ii. Implicit Midpoint Rule (SRK): Symplectic but does not preserve Lie group 
properties. 

iii. Crouch- Grossman (LGM): Lie group method but not symplectic. 

We consider the energy conservation properties of these integrators, and the 
extent to which they stay on the rotation group, as a function of step-size. 
Furthermore, we will explore how the computational cost scales as the step-size 
is varied. 



Discussion of numerical results. Notice that while one typically expects a sym- 
plectic method to exhibit good energy behavior, we find that the implicit mid- 
point method (SRK) does more poorly than the two methods that preserve the 
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Fig. 1: Comparison of Lie Group Velocity Verlet (LGVI) with Explicit Midpoint 
(RK), Imphcit Midpoint (SRK), and Crouch-Grossman (LGM). 
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Lie group structure — the Crouch-Grossman (LGM) and the Lie Group Veloc- 
ity Verlet (LGVI) methods. We find that preserving the Lie group structure 
and the symplectic structure simultaneously yields the best results in terms of 
energy preservation. 

Not surprisingly, the Lie group methods perform the best in terms of the or- 
thogonality error. The increase in orthogonality error as the step-size decreases 
is associated to the accumulation of round-off error, since smaller step-sizes 
result in more matrix multiplications. This phenomena can be mitirated by a 
more careful implementation that utilizes compensated summation |9|. 

The Lie Group Velocity Verlet method also exhibits the best computational 
efficiency, since it only requires one force evaluation per time-step, an advantage 
that it inherits from the vector space version of the Verlet method for separable 
Hamiltonian systems. 



3 Variational Integrators from Arbitrary One-Step Methods 

This section discusses how one can, given a p-th order accurate one-step method 
for an initial-value problem, and a g-th order accurate numerical quadrature 
formula, construct a variational integrator with order of accuracy min(p, q). 

A related effort to develop variational integrators from arbitrary one-step 
methods was introduced in [235, but this was developed using a nonstandard 
formulation of discrete variational mechanics on phase space and curve seg- 
ments which is quite an abstract and involved construction. In contrast, 
the approach proposed in this section is far more transparent, relies on the 
well- understood theory of shooting methods for boundary- value problems, and 



is developed in the standard setting of discrete Lagrangian mechanics [21[ . 



Outline of Approach. Notice that the characterization of the exact discrete La- 
grangian associated to Jacobi's solution is expressed in terms of the action inte- 
gral evaluated on a solution of a two-point boundary- value problem. A standard 
method of solving a boundary-value problem is to reduce it to the solution of 
an initial- value problem using the method of shooting. 

We obtain a computable approximation to the exact discrete Lagrangian ([1} 
in two stages: (i) apply a numerical quadrature formula to the action integral, 
evaluated along the exact solution of the Euler-Lagrange boundary- value prob- 
lem; (ii) replace the exact solution of the Euler-Lagrange boundary- value prob- 
lem by a converged shooting solution associated with a given one-step method. 

Shooting-based discrete Lagrangian. Given a one-step method \['h '■ TQ — > TQ, 
and a numerical quadrature formula f{x)dx ~ hJ2^=obif{x{cih)), with 
quadrature weights bi and quadrature nodes ~ cq < ci < . . . < c„_i < c„ = 1, 
we construct the shooting-based discrete Lagrangian, 

Ldiqo,qi;h)^hJ2" b^y), (10) 

where 
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Note that while we formaUy require that the endpoints are included as quadra- 
ture points, i.e., cq — 0, and c„ = 1, the associated weights bo, bn can be zero, 
so this is does not constrain the type of quadrature formula we can consider. 

Order of discrete Lagrangian. The order analysis of the shooting-based discrete 
Lagrangian depends critically on the global approximation properties of the 
shooting solution of two-point boundary-value problems. 

In general, the Euler-Lagrange equation is a second-order nonlinear differen- 
tial equation, and it is a standard result in the numerical analysis of the shooting 
method for nonlinear problems (see, for example. Theorem 2.2.2 of [lo|) that 
the approximation error in the solution of a boundary- value problem is bounded 
by the sum of two terms, the first of which is 0{h^), associated with the global 
error of the one-step method applied to the initial-value problem, and the sec- 
ond of which is related to the rate of convergence of the nonlinear solver used 
to determine the appropriate initial condition for the initial- value problem that 
would lead to the correct terminal boundary condition of the boundary-value 
problem. In essence, if the solution of the boundary-value problem is isolated 
and sufficiently regular, then a properly converged shooting method yields a 
solution of the two-point boundary- value problem that has error 0{hP), if the 
underlying one-step method has local truncation error 0{hP^^). 

This result allows us to obtain the following theorem on the order of accuracy 
of the shooting-based discrete Lagrangian, which in turn allows us to establish 
the order of accuracy of the associated variational integrator. 

Theorem 1 Given a p-th order accurate one-step method if", a q-th order ac- 
curate quadrature formula, and a Lagrangian L that is Lipschitz continuous in 
both variables, the associated shooting-based discrete Lagrangian (jlOp - (jll[) has 
order of accuracy mm{p, q). 

Proof By Theorem 2.2.2 of [13], a fully converged shooting solution (goi^'i^oi), 
associated with a one-step method ^ of order p, approximates the exact solution 
(901,^^01) of the Euler-Lagrange boundary- value problem with the following 
global error. 



qoi{c,h) ^ qoi{c,h) + 0{hP), 
VQiic^h) = iJol{c^h) + 0{hP). 



If the numerical quadrature formula is order q accurate, then 




[hYZi b^m^l{c^h) + o{hn,voi{c,h) + om)] + 




Ld{qo,qi;h) + 0{hP+^) + 0{h'^+^) 
Ld{qo,qi;h) + 0{h'^'<P''i^+^), 
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where we used the quadrature approximation error, the error estimates on the 
shooting solution, and the assumption that L is Lipschitz continuous. 

Remark 3 Notice that the numerical quadrature formula introduces an addi- 
tional factor of h which gives the requisite 0{hP^^) local error estimate for 
the discrete Lagrangian. Any consistent numerical quadrature formula would 
introduce this factor of h, since the integral is over the interval [0, /i]. 

Remark 4 It should be noted that the prolongation-collocation variational inte- 
grators introduced in flil] can be viewed as a special case of the shooting-based 
variational integrator. In particular, this involves the collocation method on the 
prolongation of the Euler-Lagrange vector field as the one-step method, and 
Euler-Maclaurin formula as the numerical quadrature method. 

More general quadrature formulas. While we have primarily discussed numeri- 
cal quadrature formulas that only depend on the integrand, one could consider 
more general quadrature formulas that depend on derivatives of the integrand, 
such as Gauss-Hermite quadrature. This would require information about the 
second- and higher-derivatives of q, which can be obtained by considering the 
prolongation of the Euler-Lagrange vector field. It is easy to show that pro- 
longations of the vector field can be used to express all higher-derivatives of q 
in terms of q, q, and derivatives of the Lagrangian L. Therefore, the shooting 
method would provide the necessary information (q, q) at the quadrature points 
to deduce the higher-derivatives {q, q^-^\ . . .), which in turn would allow the use 
of more general quadrature formulas. 

Symmetric shooting-based discrete Lagrangians. We now show that self-adjoint 
one-step methods and symmetric quadrature formulas yield self-adjoint shooting- 
based discrete Lagrangians. This guarantees that the resulting variational inte- 
grator is symmetric, and hence has even order of accuracy. 

Proposition 3 Given a self-adjoint one-step method Wh, and a symmetric 
quadrature formula (c,; -t- c„_i = 1, 6^ = bn^i), the associated shooting-based 
discrete Lagrangian is self- adjoint. 

Proof The discrete Lagrangian is given by, 

En 
. hL{q\v'-), 
2 — 

where the sequence (Q^ v*) satisfies 

The adjoint discrete Lagrangian is given by, 

En 
. 6jL(g',u'), 
2 — 

where the sequence (^^^^*) satisfies 
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Since is self- adjoint, =^h ~^-h' ^^'^ hence, 

Also, Ci + Cn-i = 1 implies that (cj+i — Cj) = (c„_i — c„_i_i). Prom these 
two properties, it is easy to see that the two sequences are related by = 
g', and = Indeed, direct substitution allows us to rewrite {q^,v'^) = 

^(c.+,-c.)h{St+\i'+^) as {q^-\v^-') = >Z'(e„_.-c„_._i)?.(9""'"''^""'"')' which 
is equivalent to = '^[ci+-,-ci)h{<l\v''). Then, 

L:iqo,qi;h) = -(-/i)^" hL{q\vl = hJ2" hL{q--\v--') 

bn-tL{q\v') =hy b,L{q\v') = Ld{qo, qi; h), 

which implies that the discrete Lagrangian is self-adjoint, where we used the 

fact that bi = b„-i. 

Group-invariant shooting-based discrete Lagrangians. By the discrete Noether's 
theorem, a G-invariant discrete Lagrangian leads to a momentum-preserving 
variational integrator. In the Galerkin discrete Lagrangian construction, this 
can be achieved by choosing a G-equivariant interpolation space. We will see 
that G-equivariant one-step methods play a similar role in the construction of 
G-invariant shooting-based discrete Lagrangians. 

Let ^:GxQ— ^Qbea group action of G on Q, and let : G x TQ TQ 
be the associated tangent-lifted action. A one-step method ^h. '■ TQ — > TQ is 
G-equivariant if it commutes with T<?, i.e., ° T^g = T^g ° for all g € G. 

Proposition 4 Given a G-equivariant one-step method '■ TQ — > TQ, and 

a G-invariant Lagrangian L : TQ — >• R, the associated shooting-based discrete 
Lagrangian is G-invariant. 

Proof For notational simplicity, we will denote the group action by 'Pg{q) = 
gq, and the tangent- lifted action by T<l>g{q,v) = {gq,gv). By definition, the 
shooting-based discrete Lagrangian is given by, 

En 
. biL{q\v'), 
1=0 

where the sequence satisfies 

{q'+\v'+')=^^,,^,_,,^h{q\v'), q° = qo, g" = Si- 

Also, 

En 
. biL{q\v'), 
1=0 

where the sequence (gSw') satisfies 

{q'+\v'+') = ^^,,^,_,,^h{q\v% f=gqo, € = gqi- 
By G-equivariance of the one-step method •Z'/^, we have, 



{gq'+\gv'+^) = 1'(c.+,-cM9q''9v'), 
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where gq^ = gqo and gq" ~ gqi. From this, we conclude that the two sequences 
are related by = {gq^,gv^). Hence, 

En ^ — 
hiL{gq\gv'') ^ b,L{q\v^) ^ La{qo,qi)- 

where we used the G-invariance of the continuous Lagrangian. 



3.1 Implementation issues 

While one can view the implicit definition of the discrete Lagrangian separately 
from the implicit discrete Euler-Lagrange equations, 

Po = -DiLdiqo,qi;h), pi = D2Ld{qo,qi;h), 

in practice, one typically considers the two sets of equations together to implic- 
itly define a one-step method: 



E7l 

. hLWy), (12a) 
1—0 

{q'+\v'+^) = ^^,^^^_,^)n{q\^l. f = 0,...n-l, (12b) 

q° = qo, (12c) 

9" = 91, (12d) 

Po = -DiLd{qo,qi; h), (12e) 

Pi^ D2Ld{qo,qi;h). (12f) 



Equation count. The unknowns in the equations are (9*,w')iLo, 90j ^ij and po, 
Pi . If the dimension of the configuration space Q is to, then there are 2{n + 3)to 
unknowns. Equations (|12cp - (|12f|) yield 4to conditions, the propagation equa- 
tions (jl2bp give 2nTO conditions, and the initial conditions {qo,po) give the last 
2to conditions. While it is possible to solve this system directly using a nonlin- 
ear root finder, like Newton's method, another possibility is to adopt a shooting 
approach, which we will now describe. 

Shooting-based implementation. Given initial conditions (goiPo), we let q" = qo, 
and guess an initial velocity Using the propagation equations (g'"*"^, v'"*"^) = 
'^(ci+i-ci)h{q^,v'^), we obtain (g',w*)"^i. Then, we let qi = q^ , and compute 
Pi ~ D2Ld{qo, qi', h). However, unless the initial velocity v'^ is chosen correctly, 
the equation pq = —DiLd{qo,qi',h) will not be satisfied, and one needs to 
compute the sensitivity of —DiLd{qo,qi]h) on w", and iterate on u° so that 
Po — —DiLd{qo,qi\h) is satisfied. This gives a one-step method (gojPo) ^ 
{qi,pi). In practice, a good initial guess for w*^ can be obtained by inverting the 
continuous Legendre transformation p = dL/dv. 
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3.2 Generalizations 

Hamiltonian Approach. It might be preferable to express the variational inte- 
grator in terms of the Hamiltonian, and an underlying one-step method for 
Hamilton's equations. We first note that, 

L{q,v) = pv-H{q,p)\^^gjj^g^. 

This leads to the following discrete Lagrangian, 

where 

= ^(,^^,_,^),M\p1, q' - 90, - qi- 

While it is possible, as before, to solve this together with the implicit discrete 
Euler-Lagrange equations as a nonlinear system, it can also be solved using a 
shooting method. 

Shooting-based implementation. Given initial conditions {qo,po), we let - go, 
guess an initial momentum p", and let v'^ — dH / dpijf ,p'^). Using the propa- 
gation equations (9'^^,^'^^) — ^{ci+i-ci)h{q^ >P^) and the continuous Legendre 
transformation = dH/dp{q^,p'^), we obtain (q*, Then, we let qi = 

g", and compute pi — D2Ld{qo, qi; h). As before, unless p° is correctly chosen, 
Pq = —DiLd{qa, qi', h) will not be satisfied, and one needs to compute the sensi- 
tivity of —DiLd{qo,qi', h) on and iterate on p° so that pa = —DiLd{qa, qi; h) 
is satisfied. This gives a one-step method {qo,po) ^ 

Type II Hamiltonian Approach. In some instances, for example, when the Hamil- 
tonian is degenerate, the Type I generating function based approach to varia- 
tional integrators using discrete Lagrangians is not applicable. A construction 
based on Type II generating functions, which corresponds to a discrete Hamil- 
tonian ll|;|l7|, can be adopted. This is based on the exact discrete Hamiltonian, 



-,E, 



h 



{qo.pi; h) = p{h)q{h) - / [p{t)v{t) - Hiq{t),p{t))l=dH/Op ^t, 



where {q(t),p{t)) is a solution of Hamilton's equations satisfying the boundary 
conditions q(0) — go, p{h) — pi. A computable discrete Hamiltonian can be 
obtained by taking 

H+{qo,pi;h) ^p^q'' -hY^^^^hlp'v' ~ H{q\p'%.^QH/dp(q\p-), 

where 

(g*+\p'+i) = >^'(c,+,-c,)/.('?^P^), q° = 90, 

The discrete Hamiltonian then yields a symplectic map via the discrete Hamil- 
ton's equations, 

qi = D2Hj{qo,pi), po = DiHj{qo,pi). 
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Shooting-based implementation. Given initial conditions {qo,Po), we let q'^ = qo, 
guess an initial momentum and let v'^ — dH/dp{q^,p'^). Using the propa- 
gation equations {q'^'^^ , p^~^^) = ^(ci+i-ci)h{q^ tP^) and the continuous Legendre 
transformation = dH/dp{q^,p^), we obtain Then, we let pi = 

p", compute qi = D2H^ {qQ,pi), and iterate on until po = DiH^{qo,pi). 



3.3 Examples 

Planar Pendulum. A planar pendulum of mass m — 1 with the massless rod of 
length / = 1 is a Hamiltonian system for which the equations of motion are 

q^p, p=-smq. 

The total energy of the system is given by the Hamiltonian H{q,p) — ^p^ — 
cos{q). 

The plots in Figure [21 represent the comparison of the standard implicit 
midpoint integrator, MID, which is known to be symplectic and the shooting- 
based variational integrator, SVIMID. Following the construction outlined in 
Section 6, we employ the midpoint rule for the initial- value problem solved by 
shooting, combined with the trapezoidal quadrature rule to obtain a computable 
discrete Lagrangian. The latter leads to a symmetric method, which is a 2nd 
order symplectic integrator. 




Fig. 2: Pendulum. Comparison of a shooting-based variational integrator 
(SVIMID) and the Imphcit Midpoint (MID). 
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Fig. 3: Pendulum. Comparison of a shooting-based variational integrator 
(SVIRK4) and the symplectic Runge-Kutta (SRK4). 



In Figure |3l we compare the performance of the shooting-based variational 
integrator constructed from the explicit Runge-Kutta method and the Simp- 
son quadrature rule with the two-stage symplectic Runge-Kutta method. Both 
methods have global error of order 4. The shooting-based method SVIRK4 ex- 



hibits a slightly larger error in energy (cf. Figure 3(c) ), but it nevertheless stays 
bounded for large time-intervals. 

Simple Harmonic Oscillator. We consider a harmonic oscillator system de- 
scribed by the equations 

q = p, P=-q- 

The total energy of the system is given by the Hamiltonian H{q, p) — ^p^ + ■ 
In Figure SI we compare two shooting-based variational integrators of order 
4. The method sviRK4Sim, used above for the planar pendulum, involves the 
4th order explicit Runge-Kutta integrator and the Simpson quadrature. The 
method sviRK4EM also employs the 4th order explicit Runge-Kutta method, 
but instead uses the Euler-Maclaurin quadrature, which involves higher-order 
derivatives of the integrant. In this case, the Euler-Maclaurin formula was trun- 
cated after the second term, and thereby requires the first-derivative of the 
Lagrangian function. 

The energy behavior for both sviRK4Sim and sviRK4EM is presented in 



Figure 4(c) We note that both the global error and the energy error is generally 



better using Simpson's rule, as compared to the Euler-Maclaurin formula, even 
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time step cpu time 

(a) Global error vs. time step. The dotted line (b) Global error vs. CPU time, 

is the reference line for the exact order. 
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(c) Energy error, step-size h = 0.2. 



Fig. 4: Simple harmonic oscillator. Comparison of the variational integrator 
sviRK4Sim based on Simpson quadrature and sviRK4EM based on Euler- 
Maclaurin quadrature. 



though the underlying one-step method is identical. This indicates that the 
specific choice of quadrature formula can have a significant effect on the global 
error and energy error properties of the constructed variational integrator. 



4 Conclusions 

We presented two general techniques for constructing discrete Lagrangians: (i) 
the Galcrkin approach, which depends on a choice of a quadrature formula, 
and a finite-dimensional function space; (ii) the shooting-based approach, which 
depends on a choice of a quadrature formula, and a one-step method. 

The order of approximation and momentum-conservation properties of a 
variational integrator are related to the order of approximation and the group- 
invariance of the discrete Lagrangian, respectively. This results in a substantial 
simplification in the analysis of variational integrators, since it is easier to verify 
the approximation and group-invariance properties of the discrete Lagrangian 
than it is to directly verify the order of accuracy and momentum-conservation 
properties of the associated variational integrator. 

For Galerkin variational integrators, the group-invariance of the discrete 
Lagrangian can further be reduced to the group-cquivariancc of the finite- 
dimensional function space. For shooting-based variational integrators, the or- 
der of the discrete Lagrangian is related to the order of the quadrature formula 
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and onc-stcp method, and the group-invariancc of the discrete Lagrangian is 
related to the group-equivariance of the one-step method. Furthermore, the 
shooting-based implementation allows the variational integrator to partially 
inherit the computational efficiencies of the underlying onc-stcp method. In 
particular, a shooting-based variational integrator constructed from an explicit 
one-step method will be more computationally efficient than one based on an 
implicit one-step method. 

These two approaches provide an explicit link between the construction of 
variational integrators, approximation theory, and one-step methods for ordi- 
nary differential equations. In particiilar, this allows one to leverage existing 
theoretical results and techniques in approximation theory and the numerical 
analysis of time-integration methods in the construction and analysis of varia- 
tional integrators. 
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